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Abstract 

A mean field approximation is employed to derive a master equation suitable for self-interacting 
baths and strong system-bath coupling. Solutions of the master equation are compared with exact 
solutions for a central spin interacting with a spin-bath. 
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I. INTRODUCTION 



The mostly condensed phase environments native to proposed technologies such as molec- 
ular electronics, laser control of chemical reactions, and quantum computing require a reex- 
amination of the problem of decoherence. The presence of intra-environmental coupling, for 
example, requires generalization beyond the uncoupled oscillator baths commonly assumed 
as a starting point in older theories^]. Environmental memory effects and the generally 
non-perturbative nature of condensed phase system-bath interactions complicate matters 
even further. 

In this manuscript we use a non-perturbative mean field approximation to derive a non- 
Markovian master equation for systems interacting with coupled baths. The master equation 
is parameter free and preserves positivity. We test the theory by direct comparison with 
exact numerical results for a model system incorporating both intra-bath coupling and strong 
system-bath coupling. Components of the theory have been published previously^, 0, Q|. 
The complete derivation and comparison with exact numerical results are presented here for 
the first time. 

Exact arguments of Nakajima and Zwanzig|j] show that master equations should be of 
linear integro-differential type (e.g., see Eq. (J2J) below). In addition, results for differential 
(i.e. Markovian) master equations suggest that dissipation should be governed by some 
suitable generalization of the generators for completely-positive-dynamical-semigroupsj^]. 
On the basis of these two conditions alone one can construct a formal master equation 
which preserves positivity and has the desired form|2j. However, this introduces an infinite 
number of unknown - and potentially time-dependent - parameters. Zwanzigp derived an 
exact integro-differential master equation but the high computational cost of calculating the 
dissipation terms precludes practical application. Hence, some level of approximation of the 



Zwanzig equation jsj], consistent with the positivity requirement, would seem a promising 
approach for development of a practical theory. Since simplifying assumptions cannot be 
made regarding the nature of the bath and since system-bath coupling may be strong, a 
mean field approximation for the system-bath interaction seems appropriate. This is the 
approach taken here. 

While the master equation derived here is uniquely adapted to systems with intra-bath 
coupling other methods not discussed could conceivably be modified for the same purpose. 
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Examples include theories which employ mixtures of quantum, semiclassical and classical 
dynamics for environmental modes ranked according to their presumed importanceQ B Q 
and approximate density functional methods [10|. Among master equation approaches only 
the Redfield theory 11] is potentially applicable, but it is valid only for weak system-bath 



coupling, violates positivity 



an incorrect long time limit 



although this can in some instances be corrected 



111]), predicts 



12j and has dissipation coefficients which diverge for finite baths. 



This last problem prevents comparison of Redfield solutions with exact data for our model 
system. Finally, uncoupled harmonic baths seem a prerequisite for stochastic wave equation 
approaches [3]. 

Section II of this manuscript outlines the derivation of the master equation with emphasis 
on the basic physical ideas. References are provided for detailed discussions of individual 
points. Section III defines a model system with strong intra-bath and system-bath coupling. 
A numerical approach for obtaining exact solutions is discussed. Section IV briefly describes 
the method employed to solve the integro-differential master equation. Finally, section V 
compares the solutions of the mean field master equation with exact solutions for the model. 

II. MEAN FIELD MASTER EQUATION 

Define a projection operator P on the total (system plus bath) density x{t) such that 

P X it) = p(t)B, (1) 

where p(t) is the system density and B is the canonical bath density. Similarly, define 
Q = 1 — P. Assuming x(0) = p(0)£>, a standard derivation then gives the Nakajima- 
Zwanzig equation p 

dp{t)/dt = -(i/h)[H,p{t)] - ! dt'K(t - t')p(t') (2) 

J o 

where H = Tt^HB} is the canonical average of the total Hamiltonian H over the states of 
the bath. The memory operator K(t) takes the form 

K{t) = Ti b {LQe~ iQLQt QLB} (3) 

where L = (l/h)[H,-] is the Liouville operator. Zwanzig did not use projection operator 



Q in his original derivation 



3, 



but his results are readily generalized and this operator has 
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been favored in the subsequent literature Q]. The operator QLQ governs the dynamics of 
system-bath interaction. 

Solutions of (J2J) can only be obtained when an explicit expression for the memory operator 
(|21) can be given. Since we cannot make simplifying assumptions about the bath or use 
perturbation theory, some sort of mean field approximation for (J3J) seems suitable. To 
approximate Q we need to understand the operator QLQ. Consider the following lemmas. 

Lemma 1: QLQ is non-Hermitian. 
Clearly QLQ is non-Hermitian if P is non-Hermitian. Consider a complete orthonormal 
basis = Ui® Vj of the Liouville-Hilbert s_pace where states U{ and Vj span the system 



and bath spaces, respectively. It then follows 



1511 that P has matrix elements 



(i,j\P\k,l) = 5 i)k Tr b {u*B}Tr b { Ul } (4) 
while the matrix elements of Pt are 

(i,j\P%l) = 6 hk Ti b {u*}Ti b { Ul B}. (5) 

Since matrix elements (jlj) and (0) differ it follows that P^ ^ P and so P is non-Hermitian. 

Lemma 2: The spectral density of QLQ is complex. 
The average spectral density of QLQ can be defined by 

Q(z) = lim -^-G 21 (z) (6) 

where G is an analytic 2x2 average Green's function 

(ri[r, 2 + (z- iA)(z* + iA^)}- 1 ) ((z - iA)[r, 2 + (** + iA^){z - iA)]' 1 ) 
((z* + iA)[rf + (z- iA)(z* + at)]" 1 ) (-r][v 2 + (z* + iA){z - iA)Y l ) 

Here rj is some real parameter and A = QLQ. [If z — x + iy then x and y denote imaginary 
and real parts of the eigenvalue of QLQ.] The angle brackets denote an average over the 
Liouville-Hilbert space i.e., for any F, 



( G n 


G l2 \ 


-( 


[g 21 


G22 J 





(F) = JmJl/mn)J2Y,( i ,j\ F \h3) (7) 
i=\ j=i 

where states denote a complete set. Defining 

G° = 




it can be shown that G satisfies the Dyson equation 

G = G° + G°m, (8) 
where £ is a self-energy which to lowest order in G is 

S = (H) + (HGH) - {H)G{H) + . . . (9) 

with 

n 



1 %A 



\ 



iA^ 



Solving (jHJ) with trucated at first order in G (self-consistent Born approximation), and 
using Eq. ©, one can show 4] that the spectral density is uniform inside an ellipse 

x 2 y 2 1 

[(AAt) - (A4)] 2 + [(AAt) + (AA)] 2 = JAM) (10) 

and zero elsewhere. [Simplified formulas for parameters (*/L4^) and (AA), suitable for 
computational use, are given in Appendix A.] 

Thus, QLQ is non-Hermitian and its spectrum is complex in general. It then follows that 
for t > we may write 

e -iQLQt = ^e-^V-is*!^)^! (11) 

3 

where u>j and jj are the real and imaginary parts of an eigenvalue of QLQ and \<pj) and ($j| 
are the associated right and left eigenvectors. Consequently, the memory operator (jSJ) can 
be written as 

K{t) = J2e-^ t e-^ t Tr b {LQ\ ( p 3 )(^\QLB}. (12) 

3 

Obviously, for a large bath a great many terms will contribute to the sum in (fT2"|) . This 
suggests the possibility of replacing K{t) by its average (in the sense of Lemma 2). 

Replacing (fT2"j) by its average, and assuming that the statistics of the eigenvalues are 
independent of the eigenvectors, gives 

K(t) = (e^e^)(£Tti,{LQ\4>i)^i\QLB}) (13) 

j 

= (cos(c^)e~ 7 *)Tr b {LQ LB} (14) 

where we have used the closure relation for the eigenvectors and the fact that for each uo 
there is a — uj to obtain the second equality. [We elsewhere call this the statistical resonance 



approximation^,!^}].] Defining the memory function W{t) = (cos(ut)e 7i ) it can be shownQ] 
that the spectral density defined in Lemma 2 gives 

W(t) = [1 - ^(Pt) 1 + \{ptf - ^(Ptf + ^(Pt) V (9t)2/8 (15) 

where 



p = [(A4t) - (A4)]/ v /(A4t) (16) 
q = [(A4 t ) + (A4)]/V(A4t) (17) 

are real parameters which depend on bath temperature. This memory function is positive, 
satisfies < W(t) < 1, and typically deviates little from gaussian form. 

Finally, assuming a Hamiltonian of the form H = H s + Hf, + J2n S^R^ where H s and 
denote system operators and and denote bath operators, this mean field type 
approximation for the memory operator yields a master equation 

d P {t)/dt = -(t/h)[H s + J2R,s^ P (t)] 

- (l/h 2 )J2C,,u fdt' W{t-t'){[ P {t')S v ,S») + [S v ,S»p{t'))}, (18) 

Hi/ "0 



where i? M = Ti h {R^B} and C^ v = Tr b {(R u — R U )(R^ — R^)B} denote canonical (i.e. B = 

e -H b /kT /Tlb{e -H b /kT^ 

averages and variances of bath operators. 
Thus, this mean field type approximation gives a master equation in which all parameters 
are known and can in principle be calculated. Moreover ()18)) can be shown to preserve 
positivity of the density matrixj^. Given that master equation ()18j) was obtained assuming 
a large bath, we should expect it to be most accurate in the thermodynamic limit. In the 
next few sections we show that sensible results are obtained even when the number of modes 
of the bath is small. 

III. SPIN SPIN-BATH MODEL 

Our model system represents two electronic states of an atomic impurity in a crystalline 
solid at low temperatures. Electric-dipole transitions from the excited electronic state are 



forbidden, but vibronic cou 
dissipation in the impurity 16] 



ing with phonons of the crystal can cause decoherence and 
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The crystal is represented by a number n s of coupled phonon modes. At low temperature 
the phonon modes can be roughly modeled as spin-1/2 modes (i.e., a)a — > a z and a^+a —>■ a x ) 
with frequencies sampled from the low energy acoustic modes of the Debye spectrum. 
We set a frequency cutoff at ud = 1.] With anharmonic phonon-phonon coupling effects 
included, but neglecting the zero point energies of the oscillators, our model Hamiltonian 
takes the form 

H = fa® + P*f> + Ao-f t °$ + Etf *?> + + ~ 9 t <W (19) 

where we arbitrarily chose ujq = .8288 as the frequency of the impurity, (3 = .01 is the coef- 
ficient of a small anharmonic correction, and Ao = 1 and A are the subsystem-environment 
and intra-environmental coupling constants. Terms one, two and four of represent the 
uncoupled modes of the subsystem (labeled 0) and environment (labeled 1 through n s ). The 
third term in ([19)1 couples the subsystem and environment, while the last term couples the 
environment with itself. The sigmas represent the Pauli matrices. In our units H = 1. Note 
that the environmental part of this Hamiltonian is non-integrable for A ^ 0. 
We calculated the reduced density matrix p(t) of the impurity via the formula 

p{t) = Pu(t) Pw(t) = E Pm Tr b {\^ m (t))(Mt)\} (20) 

\P0l{t) Poo(0 J "1=1 

where 

p m = exp{-e m //cT}/ ^ exp{-e//A;T}, (21) 

i=i 

e m and |m) are the energies and eigenvectors of the isolated environment (i.e. terms 4 and 5 
of Eq. (JEU), and kT is the temperature in units of energy. The notation Tr;,{| ip m (t))(ip m {t)\ } 
indicates a trace of the full density \i/j m (t)) (ip m {t) \ over the environmental degrees of freedom. 
The states \ip m (t)) are evolved via the Schrodinger equation from initial states 

|V>m(0)> = |1> ® \m) (22) 

under Hamiltonian The basis of eigenstates of the a z operators was used to represent all 
states. The states |0) and |1) represent down and up z-components of the spin, respectively. 
Thus, the subsystem state |1) in Eq. (|22|) means that the impurity is initially in its excited 
state. 
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Equations (|20|) and (|22|) represent an impurity in a thermal solid which is excited by a 
fast laser pulse just prior to time t = and then evolves while interacting with phono ns in 
the solid. 

The calculations reported here are for n s = 14 bath spins. The ARPACK linear algebra 
software[17] was used to calculate the lowest n eig = 20 energies and eigenvectors of the iso- 
lated environment. A temperature of kT = .02 was chosen such that no states with quantum 
number m higher than n e i g would be populated at equilibrium. The numerical solutions of 
the Schrodinger ordinary differential equations for \ip m (t)) were calculated using an eighth 
order Runge-Kutta routine Q]. Operations of the Hamiltonian on the wavevector were 
calculated via repeated application of Pauli matrix multiplication routines. For example 

O'l, • • • , ji, • • • , jn„ Wx ] 1^) = (jl, ■ ■ ■ , ji, ■ ■ ■ , 3n a ( 23 ) 

for all sets of j\ — 0, 1, I — 1, . . . , n s and where ji — 1 if ji — and ji — if j t — 1. Thus, 
an operation of g^> simply rearranges the components of \ip). States of the basis can be 
represented by integers j = j\ + ji2 + . . . + ji2 l ~ 1 + . . . + j ns 2 ?ls ~ 1 and since integers are 
represented in binary form on a computer, the mapping j — > j' = j\ + j 2 2 + . . . + jj2* _1 + 
• • •+Jn 3 2 ns_1 under cr^ can be calculated very simply using Fortran binary-operation system 
functions. Operations for crW and ay are also straightforward. 

We calculated four observables. The first is the subsystem entropy 

Sit) = -Tr{p(t)logp(t)} (24) 



1 / 1 + \/l - 4detfp(t)l 

--{logdetb(t)] + v/l-4det[p(t)]log V , = }, (25) 

2 l-Jl-4det{p{t)\ 



where det[p(t)] = pn(t)Poo(t) — Pw(t)Poi(t), which is initially zero because the reduced 
density of the subsystem is initially pure. The maximum value of this entropy is log 2 which 
corresponds to the state 

Pn(t) = l - = P w{t) (26) 
Pio(t) = = p i(t). (27) 

The entropy gives us a quantitative measure of decoherence and dissipation effects. We also 
calculated the expectations of the three components of the subsystem spin 

X(t) = Tr{af p(t)} = Pl0 (t) + p 01 (t) (28) 



Y(t) = Ti>(°V(t)} = i(p 10 (t) - poi(t)) (29) 
Zit) = Tr{af p(t)} = p u (t) - p 00 (t). (30) 

The Zit) component provides information about dissipation, while the X(t) and Y(t) com- 
ponents provide information about decoherence. 



IV. NUMERICAL SOLUTION OF MASTER EQUATION 

We recently developed a numerical technique for solving integro-differential equations I20I]. 
The accuracy of the method has been established for both generalized Langevin equations 
and master equations of type (fT%|) by comparison with exact solutions [20]. Basically the 
method works by converting integro-differential equations to ordinary differential equations. 

We implement the method as follows. Define a space-like time variable u and a smoothed 
density function 

X (t, u) = f(u) f dt> W{t -t' + u)p(t'), (31) 
Jo 

where f(u) is a damping function such that /(0) = 1. Direct substitution shows that pit) 
and x(t, u) satisfy ordinary differential equations 

dp(t)/dt = -{i/h)[H a +J2R^, P (t)\ 

- {l/h 2 ) C^{[ X (t, 0)S V , S,] + [S v , SjcQt, 0)]}, (32) 

d X (t, u)/dt = f{u)W{u)p{t) + - x (t, u) (33) 

du f{u) 

or more specifically for the spin-spin-bath model 

dp(t)/dt = -i^o® +00® ,p(t)} - 2C{ X (t,0) - a^ X (t,0)a^} (34) 
d X (t, u)/dt = e-^ 2 W{u)p{t) + + 2gu X (t, u), (35) 

where (3 = (3 + \qT> x and C = Aq(E^ — Y? x ). Here Y, x = X)&=i v^f* an d the overbar denotes a 
canonical average with respect to bath degrees of freedom. The parameters of the memory 
function were calculated using the formulas in Appendix A and the exact energies 
and eigenvectors of the bath Hamiltonian computed in Section III. The same data was 
used to calculate C and S x . Following Ref. |20j a damping function f(u) = e -9 " 2 with 



g = 11/ [(n — I) At] 2 was used. The differential equations were solved by denning a grid of 
points Uj = (n + l—j)At with j = 1, . . . ,n and I = int(.338n) where At = .1 is the time-step 
employed in the dynamics. Converged results were obtained for n = 50 grid points. We 
chose W(u) = W(\u\) for negative values of u. A discrete- variable |19j] matrix representation 
was employed to calculate the partial derivative with respect to u in Eq. (j^o]) . 

Finally, the ordinary differential equations ()34|) and (J33j) were integrated using an eighth 
order Runge-Kutta routine [l^ . 



V. RESULTS 



The exact entropy S(t) is plotted in Fig. 1(a) for intra-bath couplings A = 2 (solid curve), 
A = 4 (long-dashed) and A = 10 (short-dashed). While the results show a high degree of 
oscillation due to the relatively small number of bath degrees of freedom, there is a clear 
trend toward smaller entropy as the intra-bath coupling is increased. For A = 2 the entropy 
oscillates between zero and .5 with a mean of .25, while for A = 4 and A = 10 the mean 
values are roughly .08 and .006 respectively. Figure 1(b) shows the entropy predicted by 
the mean field master equation for the same values of the intra-bath coupling. The same 
trend toward lower entropy with higher A is observed. However, because of the mean field 
character of the theory no oscillations are observed. 

An explanation of this trend toward lower entropy for larger intra-bath coupling is pre- 
sented elsewhere jlfj]. 

Expectations of the components of the subsystem spin are plotted in Fig. 2 for A = 2 for 
exact (long-dashed) and mean field (short-dashed) calculations. For reference we also show 
the spin dynamics in the absence of system-bath coupling (solid curve). Clearly, the exact 
and mean field results show strong decoherence of comparable magnitude. After a short time 
the mean field X(t) and Y(t) become phase shifted from the exact results which also show 
evidence of noise. Amplified oscillations of similar character are observed in the exact Z(t) 
but are absent in the mean field solution. The decoherence-free solution is indistinguishable 
from the upper boundary of the figure. 

Similar calculations are shown in Fig. 3 for A = 4 and in Fig. 4 for A = 10. Decoherence 
of X(t) and Y(t) is incementally decreased in both exact and mean field solutions which show 
increasingly good agreement. Reduced dissipation in Z(t) is predicted by both calculations 
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FIG. 1: A = 2,4,10 




as A increases. 

Given the relatively small size of the present bath and the consequently oscillatory char- 
acter of the observables, and supposing that the mean field approximation will improve with 
larger baths, these results appear to support the general approach of treating the system- 
bath interaction in a mean field approximation. 

VI. SUMMARY 

The problem of predicting the dynamics of a system interacting with a condensed phase 
environment requires a reexamination of the assumptions commonly employed in theories of 
decoherence and dissipation. Specifically, intra-bath coupling cannot be neglected a priori, 
system-bath coupling may be strong, and memory effects may play a role. Simplifications 
employed in older theories such as modeling the bath as uncoupled oscillators, assuming the 
validity of perturbation theory in the system-bath coupling, and use of Markovian assump- 
tions must therefore be abandoned in general. This raises the issue of how an appropriate 
theory can be derived. 

In this manuscript we introduce a sort of mean field approximation in the system-bath 
coupling and use it to obtain an approximate master equation which preserves positivity. 
The predictions of the master equation are tested against exact results for a model system 
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FIG. 2: A = 2 
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consisting of a spin interacting with a spin-bath. In spite of the oscillatory character of some 
of the system observables good qualitative agreement is observed, raising the possibility that 
the master equation may prove quantitatively accurate for larger baths. We hope to soon 
further test the theory against exact results for a coupled oscillator bath using a recently 
developed exact method j21| for decomposing the iV-vibrational-mode time evolving density 
matrix (for pairwise interactions) into N one-dimensional stochastic density equations. 
The author gratefully acknowledges the support of the Natural Sciences and Engineering 
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FIG. 3: A = 4 




Research Council of Canada. 



VII. APPENDIX A 

Explicit formulas for (AA 1 ) and (AA) in terms of averages over finite basis sets of the 
Hilbert space are as follows: 

(AA 1 ) = — ^[2m s m 6 tr{iJ 2 }-2tr{iJ} 2 -8m s tr{iJ 2 B}-4tr{iJB} 2 + 2m s m 6 tr{if 2 B 2 } 

mjmf 
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FIG. 4: A = 10 
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+ 2m s tr{# 2 }tr 6 {B 2 } + 4m 6 tr{#B}tr{iJB 2 } + 8trt,{tr s {iJ} 2 B} + 4m s tr s {tr 6 {#B} 2 } 
+ 4m s tr s {tr 6 {#B 2 }tr 6 {#}} - 4tr{#}tr{#B 2 } - 2m 6 tr 6 {tr s {#B} 2 } 

- 2tr 6 {tr s {#} 2 }tr 6 {B 2 } - 4m s m b tr s {tr b {HB}tr b {HB 2 }} 

- 4m s tr s {tr 6 {#}tr 6 {#B}}tr 6 {B 2 } + 2m s m b tr s {tr fe {#B} 2 }tr 6 {B 2 } 

- 2m fe tr{#B} 2 tr 6 {B 2 } + 4tr{#}tr{#B}tr 6 {B 2 }] (36) 
(AA) = -^-^[2m s m b ti{H 2 }-2ti{H} 2 + AtT b {tT s {H} 2 B} + 2m s tr s {tT b {HB} 2 } 
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- 4m s tr{# 2 B}-2tr{i/B} 2 ]. (37) 

Here the traces are over the normal Hilbert space rather than the Liouville-Hilbert space. 
Specifically, the finite basis (m s basis functions) trace over subsystem degrees of freedom is 
denoted tr s {-}, the finite basis (m& basis functions) trace over reservoir degrees of freedom 
is denoted tr&{-}, and the complete trace over the finite basis (m s x m& functions) is denoted 
tr{-}. It is essential that all matrices be represented in the finite basis before the calculations 
for ()36|) and (}3*Tj) are carried out. The size of the finite basis should be chosen so that higher 
energy states are unpopulated at the given temperature. 
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